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Abstract 

Direct  internal  reforming  in  solid  oxide  fuel  cell  (SOFC)  results  in  increased  overall  efficiency  of  the  system.  Present  study  focus  on  the 
chemical  and  electrochemical  process  in  an  internally  reforming  anode  supported  SOFC  button  cell  running  on  humidified  CH4  (3%  H2O).  The 
computational  approach  employs  a  detailed  multi-step  model  for  heterogeneous  chemistry  in  the  anode,  modified  Butler- Volmer  formalism  for  the 
electrochemistry  and  Dusty  Gas  Model  (DGM)  for  the  porous  media  transport.  Two-dimensional  elliptic  model  equations  are  solved  for  a  button 
cell  configuration.  The  electrochemical  model  assumes  hydrogen  as  the  only  electrochemically  active  species.  The  predicted  cell  performances 
are  compared  with  experimental  reports.  The  results  show  that  model  predictions  are  in  good  agreement  with  experimental  observation  except  the 
open  circuit  potentials.  Furthermore,  the  steam  content  in  the  anode  feed  stream  is  found  to  have  remarkable  effect  on  the  resulting  overpotential 
losses  and  surface  coverages  of  various  species  at  the  three-phase  boundary. 

©  2006  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  possibility  of  running  hydrocarbons  (HCs)  with  or  with¬ 
out  pre-reforming  make  solid  oxide  fuel  cells  (SOFCs)  unique 
among  all  fuel  cells.  Though  the  direct  use  of  HCs  is  desirable 
from  thermodynamic  view  point,  the  breakthrough  depends  on 
the  continuous  stable  operation  of  the  cell  without  degradation. 
This  is  highly  challenging  on  conventional  Ni  based  anode  cer¬ 
mets  because  of  coking  propensity.  There  are  many  reports  on 
the  operation  of  SOFCs  running  on  H2  and  CH4  [1-5].  Works 
on  SOFCs  running  on  higher  hydrocarbons  are  also  reported 
[6-8].  Liu  and  Barnett  [5]  reported  SOFC  (LSM-YSZ/YSZ/Ni- 
YSZ)  running  on  CH4  and  natural  gas.  The  polarization  curves 
reported  show  identical  behavior  for  both  fuels.  It  could  be  im¬ 
puted  to  the  fact  that  CH4  is  the  major  constituent  of  natural 
gas.  Natural  gas  is  an  ideal  fuel  for  SOFCs  because  of  the  wide 
spread  availability  and  distribution  infrastructure. 
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The  present  work  mainly  focus  on  the  analysis  of  detailed 
chemical  processes  within  the  anode,  electrochemical  processes, 
and  losses  of  SOFCs  running  on  CH4  rich  fuel  under  internal 
reforming  conditions. 

2.  Modeling 

There  are  many  reports  in  general  on  the  modeling  and 
simulation  of  SOFC  processes.  Most  of  these  models  vary  in 
the  assumptions  made  in  deriving  the  model  equations  and 
the  dimensionality  of  the  problem.  Simple  one-dimensional 
models  to  three-dimensional  stack  simulations  have  been  re¬ 
ported  [9-15].  In  a  former  study  [9]  we  reported  the  op¬ 
eration  of  planar  SOFC  running  on  pre-reformed  CH4  fuel. 
The  model  assumed  plug  flow  in  the  fuel  and  air  channels 
and  one-dimensional  transport  in  the  electrode  structures  us¬ 
ing  Dusty  Gas  Model  (DGM).  The  electrochemistry  was  im¬ 
plemented  using  a  modified  Butler- Volmer  setting  based  on 
elementary  charge  transfer  chemistry,  and  the  anode  chem¬ 
istry  was  modeled  by  a  multi-step  heterogeneous  reaction 
mechanism. 
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Nomenclature 

As 

specific  catalytic  area  (m-1) 

(ip 

particle  diameter  (m) 

ne 

Ukj 

effective  binary  diffusion  coefficient  (m2  s-1) 

ne 

uk,  kn 

effective  Knudsen  diffusion  coefficient  (m2  s-1) 

nDGM 

ukl 

DGM  diffusion  coefficient  matrix 

E 

potential  (V),  activation  energy  (J  mol-1) 

E 

Faraday  constant  (C  mol-1) 

F 

force  (kgm-2  s-2) 

i 

current  density  (A  cm-2) 

io 

exchange  current  density  (A  cm-2) 

J 

diffusion  flux  (kg  m-2  s-1) 

Kr,Kg 

,KS  number  of  reaction,  gas-phase  species,  surface 

species 

l 

length  (m) 

P 

pressure  (Pa) 

r 

electrochemical  rate  (molm-2  s-1) 

rv 

pore  radius  (m) 

R 

gas  constant  (J  mol-1  K-1) 

T^tot 

area  specific  resistance  (Q  m2) 

s 

species  production  rate  (molm-2  s-1) 

T 

temperature  (K) 

l) 

velocity  vector  (ms-1) 

W 

molecular  weight  (kg  mol-1) 

[X] 

concentration  (molm-3) 

Yk 

mass  fraction 

Greek  symbols 

a 

permeability  (m) 

p 

symmetry  factor 

Sij 

Kronecker  delta 

0 

overpotential  (V) 

jl 

viscosity  (kgm-1  s-1) 

V 

stoichiometric  coefficient 

0 

conductivity  (S  m-1) 

X 

tortuosity 

€ 

porosity 

Subscript 

a 

anode 

c 

cathode 

e 

electrolyte 

g 

gas 

k 

species  index 

Here  we  resort  to  model  the  button  cell  experiments  reported 
by  Liu  and  Barnett  [5].  The  experiment  describes  the  opera¬ 
tion  of  an  anode  supported  button  cell  with  97  vol%  CH4  and 
3  vol%  H2O  as  the  fuel  stream  and  air  as  oxidant.  The  model 
presented  here  solves  the  velocity  v,  pressure  p  and  the  species 
mass  fractions  Yk  as  a  function  of  axial  and  radial  positions, 
and  the  electrochemical  model  parameters.  A  schematic  repre¬ 
sentation  of  button  cell  configuration  is  shown  in  Fig.  1,  and  a 


Fig.  1.  A  three-dimensional  schematic  representation  of  button  cell  configura¬ 
tion  (a)  and  the  two-dimensional  configuration  used  for  simulation  (b). 


two-dimensional  plane  of  which  is  used  as  the  model  geome¬ 
try.  The  anode  is  a  20  mm  in  diameter  and  0.5  mm  thick  porous 
membrane.  Since  the  experimental  report  does  not  give  a  detailed 
description  of  the  flow  configuration,  the  simulation  study  as¬ 
sumes  the  inlet  fuel  pipe  to  be  7  mm  in  diameter  and  ends  5  mm 
above  the  anode.  The  following  subsections  describe  various 
models  used  in  the  present  study. 


2.7.  Electrochemistry 

The  electrochemical  model  assumes  that  H2  is  the  only  elec- 
trochemically  active  species,  and  charge  transfer  reactions  occur 
only  at  the  interface  formed  by  the  electrocatalyst,  electrolyte, 
and  gas-phase  known  as  three-phase  boundary  [16].  In  the  case 
of  composite  electrodes  this  reaction  zone  can  spread  out  into 
the  electrode  usually  of  the  order  of  few  ~10  [xm  [1,17].  This 
is  a  small  fraction  of  the  total  thickness  of  the  anode  in  case 
of  an  anode  supported  cell,  and  hence  we  can  safely  assume 
that  the  charge  transfer  reaction  is  happening  at  the  well  defined 
three-phase  boundary.  However,  it  should  be  pointed  out  that 
increasing  the  three-phase  boundary  length  is  a  technology  aim 
[18]. 
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The  irreversibility  associated  with  the  cell  processes  lead  to 
various  overpotential  losses  such  as  activation,  ohmic  and  con¬ 
centration  losses  [19].  The  operating  cell  potential  is  dependent 
on  the  overpotential  losses  and  can  be  written  as  a  function  of 
local  current  density  i. 

Ecell  —  ^rev  —  ^a(0  —  l^c(OI  —  Oohm(0  —  OconciO  (1) 


and  for  oxygen  reduction 

.  _  .*  (po2/pbf/A 

'°  l°2l  +  (Po2/P*of/2 


where  and  *o2  are  adjustable  parameter  and  expressed  as  a 
function  of  temperature  by 


where  r]a  and  Oc  are  the  activation  losses  at  the  anode  and  cathode 
side,  respectively,  r]ohm  the  ohmic  overpotential,  and  r]conc  is 
the  concentration  overpotential.  Since  porous  media  transport  is 
modeled  in  detail,  the  concentration  overpotential  is  not  treated 
explicitly.  Erev  in  Eq.  ( 1 )  is  given  by  well  known  Nernst  equation. 


Erev  —  E®  + 


RT 

2F 


In 


1/2 

EH2,aEp2,c 

EH20,a 


lk  =  kn2  exp 
and 


i*o2  =  ko2  exp 


(10) 


(11) 


Expressions  for  p ^  and  p q9  can  be  found  in  Ref.  [9] . 


In  the  above  equation  E°  is  EMF  at  standard  pressure.  Eq. 
(2)  can  be  substituted  back  in  Eq.  (1).  The  ohmic  losses  in  Eq. 
(1)  is  defined  as 

Oohm  —  RioP  (3) 


where  Rto t  is  the  total  area  specific  resistance.  However,  the  resis¬ 
tance  offered  by  the  anode  and  cathode  materials  are  negligible 
compared  to  the  electrolyte  resistance  in  modem  cells.  There¬ 
fore,  in  the  following  analysis  only  the  resistance  contributed  by 
the  electrolyte  RQ  is  considered,  which  is  defined  as 

4 

*e=-,  (4) 

where  /e  is  the  thickness  of  the  electrolyte  and  cre  is  the  electrolyte 
conductivity  defined  as  a  function  of  temperature 


( re  =  3.34  x  104  exp 


10300  \ 


The  above  equation  is  estimated  for  Y SZ  electrolyte  material  and 
gives  an  ionic  conductivity  of  2.26  S/m  at  1073  K  and  0.086  S/m 
at  1273  K,  which  is  consistent  with  the  reported  values  in  Ref. 
[19].  The  functional  relation  between  the  activation  losses  and 
current  density  is  described  by  a  Butler- Volmer  formalism.  For 
hydrogen  oxidation  it  takes  the  form 


i 


=  fi) 


exp 


^(l+)3a)^ 


—  exp 


PcFjhY 

RT  J. 


and  for  oxygen  reduction 


i 


=  io 


—  exp 


PcFrjA 

RT  J 


where  io  is  the  exchange  current  density  and  /3  is  the  charge 
transfer  coefficient.  A  detailed  derivation  of  modified  Butler- 
Volmer  equations  are  given  elsewhere  [9].  In  the  present  study 
the  exchange  current  density  is  expressed  as  a  function  of  tem¬ 
perature  and  the  partial  pressures  of  the  reactants  and  products. 
The  exchange  current  density  for  hydrogen  oxidation  and  oxy¬ 
gen  reduction  is  given  by  Eqs.  (8)  and  (9),  respectively  [9]. 


.  _  .*  (/W/>h2)1/Vh2o)3/4 
10  ~lH2  1  +(PU2/P*u2)l/2 


2.2.  Heterogeneous  chemistry 


Within  a  SOFC,  anode  is  the  component  where  the  reform¬ 
ing  and  shift  reaction  proceeds.  There  is  no  significant  hetero¬ 
geneous  chemistry  proceeding  in  the  cathode  and  electrolyte.  In 
this  work,  as  reported  in  Ref.  [5]  Ni-YSZ  is  considered  to  be  the 
anode  cermet.  Ni  is  an  effective  catalyst  for  steam  reforming  of 
hydrocarbons.  However,  there  are  issues  related  to  coking  un¬ 
der  low  steam  to  carbon  (s/c)  ratio.  Kinetics  of  steam  reforming 
over  Ni/MgAl204  based  catalysts  have  been  studied  by  Xu  [20] . 
Studies  on  internal  reforming  and  related  issues  in  SOFC  anodes 
have  been  reported  by  several  others  [21-24].  Hecht  et  al.  [25] 
reported  a  multi-step  heterogeneous  reaction  mechanism  for  the 
reforming  of  CH4  on  Ni  catalysts. 

In  the  present  work  the  chemistry  is  handled  by  this  multi-step 
heterogeneous  reaction  mechanism  for  Ni  catalysts.  The  mech¬ 
anism  consists  of  42  reactions  among  6  gas-phase  species  and 
12  surface  adsorbed  species.  However,  in  the  study  reported  by 
Hecht  et  al.  [25],  the  mechanism  was  evaluated  only  at  800  °C. 
In  this  work  we  use  an  extended  version  of  the  mechanism  eval¬ 
uated  for  temperatures  between  220  and  1700  °C.  Though  the 
mechanism  is  elementary  in  nature  it  covers  the  global  aspects 
of  reforming,  water-gas  shift  and  Bouduard  reaction.  Most  of 
the  reactions  in  the  mechanism  are  expressed  in  Arrhenius  rate 
form  and  some  are  dependent  on  the  surface  coverage.  The  rate 
constants  are  expressed  as 


(12) 


where  4  is  the  rate  constant  for  the  ith  reaction,  pki  and  €ki 
the  parameters  modeling  the  coverage  dependency  of  rate  con¬ 
stant  and  0k  is  the  surface  coverage  of  kth  species.  The  rate  of 
production  of  each  specie  is  then  given  by 


Sk 


Kr  Kg+Ks 

n  wiki 


i= 1  k=  1 


(13) 


where  Kr  is  the  number  of  reactions,  Kg  and  Ks  the  num¬ 
ber  of  gas-phase  species  and  surface  species,  respectively, 
vjd  the  difference  in  stoichiometric  coefficients  of  products 
and  reactants,  vrki  the  stoichiometric  coefficients  of  products 
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and  [X]k  is  the  concentration  of  the  kth  chemical  species. 
The  complete  heterogeneous  chemistry  model  is  given  in 
Table  1.  (The  complete  mechanism  can  be  downloaded  from 
www.detchem.com/mechanisms.) 

At  the  three-phase  boundaries  Sk  also  includes  the  electro¬ 
chemical  production  rates  for  H2,  H2O  and  O2  which  are  given 
by 

rH2  =  -^,  rH2o=^,  r°2  =  -^  d4) 

2.3.  Fluid  transport 


The  flow  field  is  resolved  on  a  two-dimensional  domain  with 
axial  and  radial  coordinate  as  independent  variables  and  under 
isothermal  conditions.  The  flow  field  is  computed  by  solving  the 
equations  of  continuity 

V  •  (pv)  =  (15) 

and  momentum 

2 

V  •  ( pvv )  =  —V p - V  •  (/xVn)  +  V  •  [/jl(Vv  +  VnT)]  +  F. 

(16) 

In  Eq.  (15)  Sm  appears  only  at  the  three-phase  boundary  and  is 
the  net  production  of  all  species  given  by 

Sm  =  '52ikAsWk.  (17) 

k=l 

In  Eq.  (16)  F  is  zero  for  the  plain  media  and  in  the  porous 
media  it  is  defined  by  Darcy’s  law  as 

£  =  --«,  (18) 

a 

where  fi  is  the  viscosity  and  a  is  the  permeability.  The  species 
concentrations  are  solved  by  the  species  transport  equation  de¬ 
fined  in  the  following  form: 

V  •  ( pvYk )  =  -V  •  Ok)  +  hWkAs  (19) 


where  Yk  is  the  mass  fraction  of  kth  species  in  the  mixture,  Sk 
the  molar  production  rate  of  the  species  due  to  surface  reactions 
defined  by  Eq.  (13),  Wk  the  molecular  mass  and  As  is  the  specific 
area  available  for  surface  reactions.  The  flux  in  Eq.  (19)  is 
given  by  DGM  as  stated  in  Ref.  [9] 


E  d“GMv ra  + 


nDGM 

Ukl 


[Xi] 


ne 

UIM 


(20) 


In  Eq.  (20)  first  and  second  term  on  the  right  hand  side  represent 
the  contribution  due  to  diffusive  flux  and  viscous  flux,  respec¬ 
tively.  However,  the  viscous  flux  driven  by  pressure  gradient  is 
negligible  compared  to  diffusive  flow  in  the  porous  anode  [2] 
and  is  neglected  in  the  following  calculations.  Dk/GM  are  de¬ 
fined  as  the  matrix  of  DGM  diffusion  coefficients  and  can  be 
represented  as  a  matrix  inverse  [9] 


dDGM  =  H-\ 


(21) 


where  the  elements  of  H  matrix  are 


hkl  = 


—  +  yh 

dIm  jfkDtj 


hi  +  (hi  -  1) 


Xk_ 

ne 

ukl 


(22) 


DQkl  is  the  effective  binary  diffusion  coefficient  and  D f  kn  is  the 
effective  Knudsen  diffusion  coefficient.  The  effective  Knudsen 
diffusion  coefficient  is  given  by 


4e  Is  RT 
3r  rp\Jw~k 


(23) 


where  rp  is  the  pore  radius,  6  the  porosity  and  r  is  the  tortuosity 
of  the  porous  medium.  The  permeability  in  Eq.  (18)  is  given  by 
Kozeny-Carman  relationship 


a  = 


u  p 


72t(1  -  6)2 


(24) 


where  dv  is  the  particle  diameter. 


3.  Computational  procedure 


The  flow  field  is  solved  using  the  commercial  CFD  code  FLU¬ 
ENT  [28].  However,  the  source  terms  and  fluxes  appearing  in  Eq. 
(19)  and  the  electrochemistry  model  (Eqs.  (1),  (6),  and  (7))  are 
implemented  using  user  defined  functions  (UDF)  [27] .  During 
each  iteration  the  thermodynamic  state  variables  and  the  species 
concentrations  are  accessed  from  the  solver,  which  are  in-turn 
used  to  evaluate  the  UDF  returned  values.  Velocity  (0.04  m/s) 
inlet  boundary  conditions  and  pressure  outlet  boundary  condi¬ 
tions  are  used  for  the  calculations.  Eqs.  (1),  (6)  and  (7)  forms  a 
system  of  algebraic  equations  with  /,  da  and  r]c  as  unknowns.  At 
three-phase  boundary  for  each  computational  cell  the  residual 
form  for  the  above  variables  can  be  written  as 

F(0)  =  0  (25) 

where  the  vector  0  is  ordered  as 

&  =  [6  >?a,  >?c|.  (26) 

This  equation  system  is  solved  using  damped  Newton  itera¬ 
tions.  The  Newton  solver  normally  converges  within  three  to  four 
iterations.  However,  calculation  of  residual  requires  the  evalua¬ 
tion  of  Nernst  potential  and  the  exchange  current  densities  which 
are  dependent  on  the  partial  pressures  of  H2,  H2O  and  O2  at  the 
three-phase  boundary. 


4.  Results  and  discussion 


Numerous  calculations  have  been  carried  out  to  reproduce 
the  experimental  data  reported  by  Liu  and  Barnett  [5].  The  pa¬ 
rameters  used  for  the  calculation  are  given  in  Tables  2  and  3. 
A  comparison  of  experimentally  observed  and  simulated  po¬ 
larization  curves  is  shown  in  Fig.  2,  for  the  fuel  composition  of 
97  vol%  CH4  and  3  vol%  H2O.  In  all  cases  the  model  reasonably 
well  reproduces  the  experimental  observations.  For  all  operating 
temperatures  we  observe  a  significant  drop  in  cell  potential  at  low 
current  densities,  indicating  the  dominance  of  activation  losses. 
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Table  1 

Reaction  mechanism 
S.no. 


1. 

2. 

3. 

4. 

5. 

6. 

7. 

8. 
9. 

10. 

11. 

12. 

13. 

14. 

15. 

16. 

17. 

18. 

19. 

20. 


21. 


22. 

23. 


24. 

25. 

26. 

27. 

28. 

29. 

30. 

31. 

32. 

33. 

34. 

35. 

36. 

37. 

38. 

39. 

40. 

41. 

42. 


Reaction 

H2  +  NI(s)  +  NI(s)  -*  H(s)  +  H(s) 
02  +  NI(s)  +  NI(s)  ->  O(s)  +  O(s) 
CH4  +  NI(s)  ->  CH4(s) 

H20  +  NI(s)  ->  H20(s) 

C02  +  NI(s)  ->  C02(s) 

CO  +  NI(s)  ->  CO(s) 

H(s)  +  H(s)  ->  NI(s)  +  NI(s)  +  H2 
O(s)  +  O(s)  ->  NI(s)  +  NI(s)  +  02 
CH4(s)  -*  CH4  +  NI(s) 

H20(s)  ->  H20  +  NI(s) 

C02(s)  -*  C02  +  NI(s) 

CO(s)  ->  CO  +  NI(s) 

^CO(s) 

H(s)  +  O(s)  ->  OH(s)  +  NI(s) 
OH(s)  +  NI(s)  -*  H(s)  +  O(s) 

H(s)  +  OH(s)  ->  H20(s)  +  NI(s) 
H20(s)  +  NI(s)  ->  H(s)  +  OH(s) 
OH(s)  +  OH(s)  ->  H20(s)  +  O(s) 
H20(s)  +  O(s)  ->  OH(s)  +  OH(s) 
C(s)  +  O(s)  ->  CO(s)  +  NI(s) 
CO(s)  +  NI(s)  ->  C(s)  +  O(s) 

^CO(s) 

CO(s)  +  O(s)  ->  C02(s)  +  NI(s) 

^CO(s) 

C02(s)  +  NI(s)  ->  CO(s)  +  O(s) 
HCO(s)  +  NI(s)  ->  CO(s)  +  H(s) 

^CO(s) 

CO(s)  +  H(s)  ->  HCO(s)  +  NI(s) 
HCO(s)  +  NI(s)  -*  CH(s)  +  O(s) 
CH(s)  +  O(s)  ->  HCO(s)  +  NI(s) 
CH4(s)  +  NI(s)  ->  CH3(s)  +  H(s) 
CH3(s)  +  H(s)  ->  CH4(s)  +  NI(s) 
CH3(s)  +  NI(s)  ->  CH2(s)  +  H(s) 
CH2(s)  +  H(s)  ->  CH3(s)  +  NI(s) 
CH2(s)  +  NI(s)  ->  CH(s)  +  H(s) 
CH(s)  +  H(s)  ->  CH2(s)  +  NI(s) 
CH4(s)  +  NI(s)  ->  C(s)  +  H(s) 

C(s)  +  H(s)  ->  CH(s)  +  NI(s) 
CH4(s)  +  O(s)  -*  CH3(s)  +  OH(s) 
CH3(s)  +  OH(s)  ->  CH4(s)  +  O(s) 
CH3(s)  +  O(s)  ->  CH2(s)  +  OH(s) 
CH2(s)  +  OH(s)  ->  CH3(s)  +  O(s) 
CH2(s)  +  O(s)  ->  CH(s)  +  OH(s) 
CH(s)  +  OH(s)  ->  CH2(s)  +  O(s) 
CH(s)  +  O(s)  ->  C(s)  +  OH(s) 

C(s)  +  OH(s)  ->  CH(s)  +  O(s) 


Aa  (cm,  mol,  s) 

P 

Eaa  (kJ  mol 

0.010  x  10-°0b 

0.0 

0.0 

0.010  x  10-°0b 

0.0 

0.0 

8.000  x  10-°3b 

0.0 

0.0 

0.100  X  Kr°°b 

0.0 

0.0 

1.000  X  10-°5b 

0.0 

0.0 

5.000  X  10“01b 

0.0 

0.0 

2.545  x  10+19 

0.0 

81.21 

4.283  x  10+23 

0.0 

474.95 

8.705  x  10+15 

0.0 

37.55 

3.732  x  10+12 

0.0 

60.79 

6.447  x  10+07 

0.0 

25.98 

3.563  x  10+11 

0.0 

111.27 

— 50.0C 

5.000  x  10+22 

0.0 

97.9 

1.781  x  10+21 

0.0 

36.09 

3.000  x  10+20 

0.0 

42.7 

2.271  x  10+21 

0.0 

91.76 

3.000  x  10+21 

0.0 

100.0 

6.373  x  10+23 

0.0 

210.86 

5.200  x  10+23 

0.0 

148.1 

1.354  x  10+22 

-3.0 

116.12 

— 50.0C 

2.000  x  10+19 

0.0 

123.6 

— 50.0C 

4.653  x  10+23 

-1.0 

89.32 

3.700  x  10+21 

0.0 

0.0 

50. 0C 

4.019  x  10+2° 

-1.0 

132.23 

3.700  x  10+24 

-3.0 

95.8 

4.604  x  10+20 

0.0 

109.97 

3.700  x  10+21 

0.0 

57.7 

6.034  x  10+21 

0.0 

61.58 

3.700  x  10+24 

0.0 

100.0 

1.293  x  10+22 

0.0 

55.33 

3.700  x  10+24 

0.0 

97.1 

4.089  x  10+24 

0.0 

79.18 

3.700  x  10+21 

0.0 

18.8 

4.562  x  10+22 

0.0 

161.11 

1.700  x  10+24 

0.0 

88.3 

9.876  x  10+22 

0.0 

30.37 

3.700  x  10+24 

0.0 

130.1 

4.607  x  10+21 

0.0 

23.62 

3.700  x  10+24 

0.0 

126.8 

1.457  x  10+23 

0.0 

47.07 

3.700  x  10+21 

0.0 

48.1 

1.625  x  10+21 

0.0 

128.61 

Total  surface  site  density  isA.  =  2.6xlO  9  mol/ cm2. 

a  Arrhenius  parameters  for  the  rate  constant  written  in  the  form:  k  =  AT^  exp( 
b  Sticking  coefficient. 
c  Coverage  dependent  activation  energy. 

However,  the  major  discrepancy  between  the  model  predicted 
results  and  the  experimentally  observed  data  lies  in  the  open  cir¬ 
cuit  voltages  (OCV).  The  experiment  reports  maximum  OCV  of 
1 . 17  at  800  °C.  However,  the  model  predicts  a  much  higher  OCV 
(1.55  V)  at  800  °C.  To  validate  the  model  predictions,  OCVs  are 
calculated  based  on  equilibrium  predictions  (with  and  without 
surface  carbon)  assuming  H2  oxidation  mechanism.  Experimen¬ 
tally  observed  OCVs  are  compared  with  those  predicted  by  the 
model  and  equilibrium  calculations  (Fig.  3).  In  general  though 


E/RT). 


slightly  higher,  equilibrium  composition  with  surface  carbon 
yields  OCVs  which  are  in  reasonable  agreement  with  experi¬ 
mental  observations.  The  experimental  observation  of  OCV  is 
lower  presumably  due  to  slight  gas  leakage.  However,  equilib¬ 
rium  composition  without  accounting  surface  carbon  results  in 
much  higher  OCVs.  At  higher  temperatures  model  predictions 
are  close  to  those  of  equilibrium  predictions  without  surface 
carbon.  While  at  low  temperatures  model  predictions  are  close 
to  equilibrium  calculations  with  surface  carbon.  This  leads  to 
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Fig.  2.  Voltage  and  power  density  against  current  density.  Comparison  with 
experimental  data  [5]. 


the  conclusion  that,  though  the  surface  chemistry  model  is  ca¬ 
pable  of  predicting  surface  carbon  formation,  refinement  of  the 
model  for  carbon  formation  is  required  for  accurate  predictions 
especially  at  higher  temperatures.  Nevertheless,  it  should  be  no¬ 
ticed  that  under  short-circuit  conditions  the  model  is  in  good 
agreement  with  experiments.  Liu  and  Barnett  [5]  have  analyzed 
a  number  of  possible  electrochemical  oxidation  reactions  and 
concluded  that  the  OCVs  for  the  partial  oxidation  of  C  are  in 
good  agreement  with  experimental  observations.  However,  we 
believe  that  when  enough  H2O  is  present  in  the  feed,  the  internal 
reforming  can  lead  to  H2  production  within  the  anode  and  under 
such  conditions  H2  oxidation  will  be  the  dominant  mechanism 
of  electrochemical  charge  transfer  compared  to  other  possible 
pathways. 

Figs.  4  and  5  show  the  cathodic  (r]c )  and  anodic  (77a)  activa¬ 
tion  losses  versus  current  density.  The  cathodic  activation  losses 
show  expected  behavior.  The  losses  increase  with  decreasing 
temperature  and  increasing  current  density.  However,  the  anodic 
activation  loss  shows  a  marked  difference  from  the  expected  be¬ 
havior.  The  anodic  losses  are  distinctly  different  for  high  and 
low  temperatures.  At  high  temperatures  the  losses  are  high  at 


Fig.  4.  Cathode  overpotentials  as  a  function  of  current  density  for  different 
operating  temperatures. 


low  current  densities.  At  this  point  it  is  worth  to  analyze  the 
functional  dependency  of  exchange  current  density  z'o-  The  ex¬ 
change  current  density  is  the  current  density  of  charge  transfer 
reaction  at  equilibrium  electric  potential  difference  between  the 
electrode  and  the  electrolyte  phases,  and  is  usually  a  strong  func¬ 
tion  of  species  composition  and  temperature.  A  high  exchange 
current  density  causes  the  electrochemical  charge  transfer  re¬ 
action  to  proceed  rapidly  upon  varying  the  potential  difference 
from  its  equilibrium  value.  A  careful  analysis  of  Eq.  (8)  reveals 
that  a  zero  H2O  partial  pressure  leads  to  zero  exchange  current 
density.  At  low  current  densities,  nearly  all  H2O  produced  at  the 
three-phase  boundary  (TPB)  by  electrochemical  charge  trans¬ 
fer  reaction  is  consumed  by  the  reforming  chemistry,  leading  to 
very  low  exchange  current  density  for  H2  oxidation.  This  basi¬ 
cally  requires  high  activation  overpotential  to  drive  the  electro¬ 
chemical  charge  transfer  reactions.  It  should  be  noticed  that  the 
activation  overpotential  is  the  potential  difference  above  the 
equilibrium  electric  potential  difference  between  the  electrode 
and  electrolyte  phases. 


Fig.  3.  A  comparison  of  experimentally  observed  OCVs  and  those  predicted 
by  the  model  and  equilibrium  calculations  with  and  without  surface  carbon 
formation. 


Fig.  5.  Anode  overpotentials  as  a  function  of  current  density  for  different  oper¬ 
ating  temperatures. 
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Table  2 


SOFC  parameters 


Parameters 

Values 

Anode 

Thickness  (/a,  mm) 

0.50 

Average  pore  radius  (rp,  p,m) 

0.50 

Average  particle  diameter  (dp,  p,m) 

2.50 

Specific  area  (As,  cm-1) 

1025 

Porosity  (6) 

0.35 

Tortuosity  (t) 

3.80 

Charge  transfer  coefficient  0 % ) 

0.50 

Electrolyte 

Thickness  (/e,  pan) 

25.0 

Cathode 

Thickness  (/c,  pan) 

30.0 

Average  pore  radius  (rp,  pan) 

0.50 

Average  particle  diameter  ( dp ,  pan) 

2.50 

Porosity  (6) 

0.35 

Tortuosity  (t) 

3.80 

Charge  transfer  coefficient  (/ % ) 

0.5 

0.040 -| 
0.035- 
0.030- 

2  0.025- 

O 

%  0.020- 

o> 

CO 
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Fig.  8.  Carbon  formation  as  a  function  of  current  density  for  the  operating  tem¬ 
perature  of  1023  K. 


Fig.  6.  Exchange  current  density  z'o  and  i/io  for  the  anodic  branch  as  a  function 
of  current  density  for  the  operating  temperature  of  1073  K. 


Fig.  7.  Surface  coverage  of  carbon  at  OCVs  for  different  operating  temperatures. 


Fig.  9.  Fraction  of  Ni  vacancies  and  surface  coverage  CO  at  the  three-phase 
boundary  as  a  function  of  current  density  for  different  operating  temperatures. 


Fig.  10.  Surface  coverage  of  hydrogen  at  the  three-phase  boundary  as  a  function 
of  current  density  for  different  operating  temperatures. 
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Fig.  11.  Surface  coverage  of  oxygen  at  the  three-phase  boundary  as  a  function 
of  current  density  for  different  operating  temperatures. 

However,  the  trend  is  different  at  low  temperatures,  where 
there  is  still  enough  H2O  at  the  TPB  due  to  the  low  rate  of  re¬ 
forming  and  hence  leading  to  the  expected  behavior  in  loss  po¬ 
tential.  At  high  temperatures  and  high  current  densities  plenty 
of  H2O  available  at  the  TPB  and  hence  the  anodic  losses  shows 
the  expected  behavior  at  high  current  densities.  Though  one  can 

Table  3 


Parameters  for  exchange  current  density 


Parameters 

Value 

H2  oxidation  (7^  ) 

&H9  (A  cm  2) 

2.07  x  105 

(kJmol-1) 

87.8 

O2  reduction  ( i q  ) 

&O9  (A  cm  2) 

5.19  x  104 

Eq0  (kJmol-1) 

88.6 

argue  this  behavior  is  a  limiting  case  of  the  exchange  current 
density  function  (Eq.  (8)),  it  is  worth  remembering  that  any  phys¬ 
ically  realistic  functional  formulation  of  exchange  current  den¬ 
sity  should  be  dependent  on  the  concentration  of  the  reactants 
and  products  participating  in  the  charge  transfer  chemistry  [10]. 
Since  H2  is  undoubtedly  an  electrochemically  active  species  any 
formulation  of  exchange  current  density  will  be  dependent  on 
H2O  partial  pressure  and  will  lead  to  same  behavior  at  low  cur¬ 
rent  densities  for  any  fuel  with  very  low  H2O  content.  A  plot 
of  exchange  current  density  z‘o  and  i/io  for  the  anodic  branch 
is  shown  in  Fig.  6.  It  is  quite  clear  that  the  anodic  overpotential 
follows  the  trend  of  i/io. 

Fig.  7  shows  the  carbon  deposition  at  OCVs  for  various  op¬ 
erating  temperatures.  Carbon  deposition  was  maximum  for  the 
highest  operating  temperature  of  800  °C.  Fig.  8  shows  carbon 
formation  as  a  function  of  current  density.  It  is  evident  from  the 
figure  that  the  flow  of  current  mitigates  coking.  It  is  mainly  be¬ 
cause  of  the  fact  that,  as  current  starts  to  flow  more  and  more 
H2O  is  formed  at  the  three-phase  boundary  and  hence  reducing 
the  possibility  of  CH4  cracking  on  Ni  surfaces.  This  observation 
of  high  C  deposition  at  OCV  is  in  good  agreement  with  the  ex¬ 
perimental  report  [5].  Surface  coverages  of  other  species  at  the 
three-phase  boundary  are  shown  in  Figs.  9-11.  Fig.  9  shows  the 
surface  coverage  of  CO  and  free  Ni  surfaces.  It  can  be  seen  that 
at  open  circuits  the  Ni  surfaces  are  relatively  open.  As  current 
starts  flowing  the  free  Ni  surfaces  are  mostly  covered  by  CO  and 
hydrogen  (Fig.  10).  However,  the  relative  coverage  of  hydrogen 
is  less  compared  to  CO.  The  surface  coverages  of  hydrogen  can 
result  from  the  dissociative  adsorption  of  H2O  and  CH4.  For  all 
temperature  regimes  major  species  on  the  surface  were  found  to 
be  hydrogen  and  CO.  However,  hydrogen  coverage  was  ~45% 
of  CO  coverage  for  all  cases.  Fig.  1 1  shows  surface  coverages  of 
oxygen.  It  should  be  noticed  that  oxygen  on  the  surface  results 
from  the  dissociative  adsorption  of  water.  For  all  the  species, 


Fig.  12.  Composition  of  anode  stream  at  OCVs  as  a  function  of  temperature. 
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Fig.  13.  H2  mole  fractions  within  the  anode  compartment. 


Fig.  15.  Velocity  (m/s)  profiles  within  the  anode  compartment. 


coverages  were  found  to  decrease  with  increasing  temperature. 
It  is  obvious  that  the  desorption  rates  are  higher  at  higher  temper¬ 
ature  thereby  leaving  the  Ni  surfaces  open  at  high  temperatures. 
We  believe  that  carbon  formation  and  its  subsequent  reaction 
with  H2O  (globally  stated  as  C  +  H2O  CO  +  H2)  plays  a  key 
role  at  low  current  densities,  where  CO  and  H  surface  coverages 
keep  increasing.  For  example,  in  the  case  of  750  °C  a  compari¬ 
son  of  surface  coverages  of  H  (Fig.  10),  CO  (Fig.  9),  and  C  reveal 
that  coking  is  suppressed  at  current  density  of  ^0.3  A/cm2  (Fig. 
8),  and  we  see  a  transition  in  the  trend  of  H  and  CO  coverages 
at  the  same  current  density. 


Radial  position  (m) 

Fig.  14.  H2O  mole  fractions  within  the  anode  compartment. 


The  composition  of  anode  stream  at  OCV  for  various  op¬ 
erating  temperatures  is  shown  in  Fig.  12.  Increasing  tempera¬ 
ture  is  found  to  increase  all  reaction  products.  CO2  is  omitted 
from  the  figure  due  to  its  very  low  concentration  at  OCVs.  Con¬ 
tour  plots  of  H2  and  H2O  mole  fractions  are  shown  in  Figs.  13 
and  14,  respectively.  It  is  quite  clear  that  H2  is  produced  within 
the  anode  and  H2O  is  produced  by  the  electrochemical  charge 
transfer  reaction  at  the  three-phase  boundary.  It  can  also  be  seen 


Fig.  16.  Comparison  of  surface  coverages  of  various  surface  adsorbed  species 
for  two  different  fuel  compositions  with  varying  H2O  content. 
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Fig.  17.  Comparison  of  overpotentials  for  two  different  inlet  fuel  compositions 
with  varying  H2O  content. 


that  the  very  low  flow  rate  causes  the  back  diffusion  of  H2  pro¬ 
duced  within  the  anode  into  the  fuel  inlet  channel.  For  the  case 
of  800  °C  velocity  profiles  within  the  anode  compartment  are 
shown  in  Fig.  15. 

To  understand  the  role  of  H2O  simulation  are  carried  out  with 
highly  humidified  CH4  (40  vol%  H2O).  A  comparison  of  surface 
coverages  of  H  and  CO,  and  empty  Ni  surface  is  shown  in  Fig.  16. 
Fuel  gas  with  40%  H2O  resulted  in  higher  H  coverage  and  lower 
CO  coverage  compared  to  fuel  gas  with  3%  H2O.  Fig.  17  shows 
a  comparison  of  activation  losses  for  anode  and  cathode  for  the 
two  different  feed  gas  compositions.  For  the  case  of  fuel  gas  with 
40%  H2O  activation  overpotentials  show  a  linear  behavior  with 
increasing  current  density  and  result  in  a  much  low  activation 
losses  compared  to  fuel  gas  with  3%  H2O. 

In  a  recent  work  Lin  et  al.  [26]  reported  the  exit  fuel  compo¬ 
sitions  from  a  button  cell  working  under  the  same  fuel  composi¬ 
tions  as  in  Ref.  [5].  They  found  all  the  products  to  be  increasing 
with  increasing  current  density,  though  the  equilibrium  calcula¬ 
tion  predicts  different  trend.  However,  data  is  reported  only  for 
lower  current  densities  (up  to  0.8  A/cm2).  Our  findings  also  pre¬ 
dict  an  increasing  product  composition  at  low  current  densities 
(up  to  0.5  A/cm2)  and  then  a  decreasing  H2  and  CO  concentra- 


Fig.  18.  Exit  gas  composition  as  a  function  of  current  density  at  800°  for  the 
case  of  97%  CH4  and  3%  H20. 


tions  for  further  increase  in  current  densities  (Fig.  1 8).  However, 
our  calculations  predict  more  H2O  than  CO2,  while  the  exper¬ 
iment  predicts  higher  CO2  than  H2O.  Our  prediction  of  higher 
H2O  than  CO2  is  consistent  with  the  equilibrium  predictions. 

5.  Conclusions 

In  the  present  work  a  detailed  CFD  study  of  the  chemical  and 
electrochemical  processes  in  an  internally  reforming  anode  sup¬ 
ported  SOFC  button  cell  was  carried  out.  We  have  implemented 
detailed  models  for  chemistry,  electrochemistry  and  porous  me¬ 
dia  transport  into  the  commercial  CFD  code  FLUENT.  Simula¬ 
tion  results  were  compared  with  experimentally  reported  data. 
The  comparisons  lead  to  the  conclusion  that  precise  calculation 
of  surface  carbon  formation  is  critical  for  the  accurate  predic¬ 
tion  of  OCVs  for  hydrocarbon  fuels  with  very  low  H2O  content. 
Anodic  overpotentials  showed  remarkable  difference  from  ex¬ 
pected  behavior.  Though  there  are  no  experimental  findings  to 
validate  our  findings  on  the  anodic  loss  potentials,  any  functional 
form  of  exchange  current  density  cited  so  far  in  the  literature  (e.g. 
[10,32])  would  lead  to  the  same  trend  as  reported  here.  Within 
the  framework  of  simulation,  the  well  understood  aspects  are 
the  flow  field  models  and  the  reforming  chemistry.  The  hetero¬ 
geneous  chemistry  model  used  in  this  work  is  well  validated  for 
Ni-YSZ  cermets  by  Hecht  et  al.  [25].  There  is  no  question  on 
the  validity  of  Navier  Stokes  equations  for  fluid  transport.  DGM 
have  been  used  by  many  researchers  for  porous  media  transport 
[9,29,30].  However,  electrochemistry  remains  as  one  of  the  least 
understood  aspect  of  fuel  cell  modeling.  As  done  in  this  work 
most  the  modeling  effort  uses  Butler- Volmer  formalism,  except 
some  state  space  modeling  works  [31,33].  Any  deviation  from 
physically  realistic  behavior  can  raise  questions  on  the  validity 
of  existing  electrochemical  models,  which  can  be  well  applied 
for  any  fuel  composition.  The  model  presented  here  is  well  ap¬ 
plicable  for  any  fuel  cell  configuration  and  can  be  applied  to 
understand  the  underlying  chemical  and  physical  processes  and 
hence  to  choose  the  best  operating  conditions  for  SOFCs. 
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